Lippmann-Schwinger description of multiphoton ionization 

I. A. Ivanov* 1 " and A. S. Kheifets 
Research School of Physical Sciences and Engineering, 
The Australian National University, Canberra ACT 0200, Australia 

(Dated: February 2, 2008) 

Abstract 

We outline a formalism and develop a computational procedure to treat the process of multi- 
photon ionization (MPI) of atomic targets in strong laser fields. We treat the MPI process nonper- 
turbatively as a decay phenomenon by solving a coupled set of the integral Lippmann-Schwinger 
equations. As basic building blocks of the theory we use a complete set of field-free atomic states, 
discrete and continuous. This approach should enable us to provide both the total and differential 
cross-sections of MPI of atoms with one or two electrons. As an illustration, we apply the proposed 
procedure to a simple model of MPI from a square well potential and to the hydrogen atom. 
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I. INTRODUCTION 



In recent years, the process of multiphoton ionization (MPI) of atomic and molecular 



species has been a su 
by Protopapas et al. 



2ject of intensive experimental and theoretical studies (see reviews 



l|, Lambropoulos et al. |2j, Chu and Telnov |^| and Posthumus J] 
and references therein). Rapid progress in this field has been largely driven by advance- 
ment in high-power short-pulse laser techniques. The laser intensities which may go beyond 
10 13 Wcm" 2 make it possible to observe many striking phenomena such as MPI and above- 
threshold ionization. 

Accurate theoretical description of ionization processes occurring in laser fields of such 
intensities should necessarily go beyond a simple perturbative picture. The first nonpertur- 
bative theory of MPI was proposed by Keldysh Q, Faisal Q and Reiss [7|. Their theory 
(known as KFR) treated the process of MPI as a transition of an electron from an initial 
bound state into a final state described by the classical Volkov wave function. The KFR 
approach provided simple analytical formulas for the MPI rate which were found in a qual- 
itative agreement with experiment. Various modifications of the KFR theory were made, in 
particular those accounting for the rescattering process si, j|. 

The KFR theory treated the laser field purely classically. The MPI problem can also 
be formulated in an entirely quantum form. The properties of the scattering matrix in 
ihis formulation of the MPI process have been studied starting from the works of Mower 
lfli 11 1, Gontier et al. 12, Faisal and Moloney Jackson and Swain [l5 |. 



recently, Guo and Aberg 



Jackson and Swain 115 1. More 

n 

161 ] and Guo, Aberg and Crasemann |17( put forward an MPI theory 
(referred hereafter as GAC) which treats the photoionization process as a QED scattering 
phenomenon. The emphasis in this theory was placed on a proper QED description of an 
electron interacting with the laser field (the quantum version of the Volkov states). Further 
development of the QED picture of the MPI phenomenon was made by Gao et al. Q] and 
Chen et al. [l^f who refined the original GAC theory by including the non-laser modes of 
the electromagnetic field. 

This fully QED approach, although solving the problem in principle, was found to be 
rather difficult to implement in practice even for simplest atomic targets such as one- or two- 
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electron atoms. For instance, Chen et al. [19j had to make further drastic approximations 
in order to carry out their calculation of MPI on atomic hydrogen. 

One possible solution to this problem is to use some suitable square-integrable basis to 
represent the Green functions occuring in the theory. Such approach, using ideas of the 



201 ]. Another solution is to 



complex rotation method, has been proposed by Maquet et al. 
use the field-free atomic states as building blocks of the theory. Such a choice is particularly 
advantageous when theory is applied to complex atomic systems with more than one target 
electron. Such approach within the context of the complex rotation method has been used 



22|. Description of the 



by Mercouris and Nicolaides [21[ and Nicolaides and Mercouris 
quantum mechanical evolution of an atom coupled to a laser field is also possible in the 
framework of this method j^ . 

In the present paper, we outline a quantum formalism for MPI which we intend to use 
for practical computations on complex atomic systems. The formalism is not entirely new 
and is based on the ideas expressed earlier in the literature. The main emphasis of this 
work is on the pracical implementation of the formlism and turning it into an efficient 
computational procedure. In this development we are inspired by a series of works by Burke 
and collaborators who combined the Floquet description of the laser field with the i?-matrix 



scattering theory. Following the seminal work 



24j , this approach has been success fully 



implemented for calculating the total MPI rate and the level shift in atomic hydrogen 
helium j^j], the negative hydrogen ion 22] and molecular hydrogen [2^. Most recently the 
i?-matrix Floquet theory was combined with the basis spline technique to describe the two- 

n n 

electron MPI from the helium atom in the ground [29] and excited states [30]. In addition to 
the total MPI rate, some differential cross-sections can be also calculated within the Floquet 
formalism as was demonstrated for atomic hydrogen by Potvliege and Shakeshaft (31I and 
Potvliege j3] • It is the most detailed fully differential cross-sections that are of particular 
interest to experimentalists and that we intend to evaluate in our approach. 

In the present paper we employ the operator formalism due to Goldberger and Watson 
33I ] . In this formalism, the MPI process is treated as a decay phenomenon. The partial 
decay rates and the energy level shifts are evaluated via the matrix elements of the transition 
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operator which are found by solving a coupled set of the integral Lippmann-Schwinger 
equations. In this approach the matrix elements of the transition operator should be taken 
between the field-free atomic states accompanied by an integer number of the laser photons. 

For one-electron targets, evaluation of the field-free states is trivial. For two-electron 
targets, an accurate set of target states, both discrete and continuous, can be generated by 
the so-called convergent close coupling (CCC) method. This method has been extensively 
tested for pro cesses with two electrons in the continuum such as electron scattering on atomic 

We intend to use the same 



hydrogen 3^ and low-field double ionization of helium 35, 3(| 



set of target states for MPI of He in the non-perturbative strong-field regime. 

In our approach we ignore all the processes of spontaneous emission of photons, which is 
ustified as we are interested in processes induced by strong fields. Since the work of Shirley 
37| , it is known that if one neglects spontaneous processes in the quantum description of 
the interaction of laser light and atom, then the quantum picture and the Floquet method 
of solving the time-dependent Schrodinger equation become equiavalent. Therefore, the 



24j . for example) are 



present approach and the ones based on the Floquet anzats (such as 
completely equivalent. However, recasting the MPI description in the form presented below 
may have some practical advantages, especially if we have good means of representing the 
field-free atomic states for complex atomic systems. This includes, of course, the states 
belonging to the continious spectrum. Such means do exist. As mentioned above, the CCC 
method turned out to be very efficient at representing field-free states of atomic systems. 
It was the desire to exploit fully the possibilities provided by CCC technique that largely 
motivated us to give a formulation of the MPI process presented below. 

The rest of the paper is organized as follows. In Section |H] we give a formulation of the 
MPI theory in terms of the field-free atomic states. In Section III Al we consider a model 
square- well problem, and in Section III Bl the MPI of hydrogen. We conclude in Section 
by outlining a set of problems we intend to consider in the immediate future. 
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II. NONPERTURBATIVE FRAMEWORK FOR THE DESCRIPTION OF THE 
MPI PROCESS. 



Let us consider a system which consists of a number of photons with a given frequency 
lo and momentum vector k corresponding to an incident plane-wave, and a target (atom or 
ion). We shall describe the field fully quantum- mechanically and write the Hamiltonian of 

the system as 

H = H atom + iJfidd + H int . (1) 
Here H atom and ifgeid have the usual meaning of the Hamiltonians of the atom and the field: 

N p 2 N Z N 1 

#atom = X! ^2 ~ 

i=i * j=i 't ij=i,i>j 'y 

i?flel d = iVcU (2) 

The atomic Hamiltonian is taken in a non-relativistic form. The number operator N refers 
to the laser photons only. 

The corresponding states of the system consisting of the non-interacting atom and the 
field are denoted as \a) = \a,n), where a set of quantum numbers a defines a state of 
the atom and n is the number of the laser photons. The following notations will be kept 
throughout the paper: Greek letters will be used to designate the states of a whole system 
"the atom plus external field", while the Latin letters will be used for the atomic states. 
The atomic system of units is in use with e = m — % — 1. 

The part of the Hamiltonian H int which describes the interaction of the atom and the 
linearly polarized laser field characterized by angular frequency oo, wave-vector k and polar- 
ization vector e can be written as (see e.g. Sobelman |38]) 

1 N ( . A 2 \ 

where p is the momentum operator and A is a quantized vector potential normalized to the 
unit volume: 

hire 2 

A = J e(ate lkr + a k e' ikr ), (4) 

V uj 
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Here a + and a have the usual meaning of the operators creating and annihilating a photon. 
As mentioned above, we are interested in sufficiently strong field intencities, when all pro- 
cesses of spontaneous emission of photons can be neglected. For this reason, only one mode 
corresponding to the wavevector of the incident laser filed is kept in Eq. (@J). Therefore, here 
and below, A describes the laser photons only. We also restrict ourselves with the dipole 
approximation in which the operator A does not act on the atomic coordinates. 

The matrix elements of the vector potential operator taken between the noninteracting 
states of the system "atom plus laser photons" are given by the well-known formulas (see 



e.g. Sobelman 



38|): 



~ / 2<7TC 2 Tl 

(a, n\ A • a\b, n — 1) = \ (a\e ■ a\b) 



i i 2 - 17 t \ /27r(n + l)c 2 ,_, 

(a, n\A- a\b, n + 1) = \ — — (a\e-a\b), (5) 

V 

where a is an arbirtary vector operator acting only on the coordinates of the atomic sub- 
system. 

In strong fields n ~ n — 1 ^> 1 and the coefficients in (j3J can be simplified to 



27mc 2 2n(n + l)c 2 Fc 

V— K V u "a? (6) 

where F is the electric field strength related to the energy density as F 2 = 8imuj. This leads 
to the following formulas for the matrix elements of the operator H in t'- 

F 

{a,n\H int \b,n±l} =- —(a\e-p\b), (7) 

F 2 

(a,n\H int \b,n±2) = — -(a \b) 

OUJ z 

F 2 

(a,n\H int \b,n) = ^( a \ b ) 
We shall treat the MPI process as a decay phenomenon within the framework of the quantum 



decay theory as described by Goldberger and Watson [33|. We shall be interested in the 
following process. At the moment t = the system "atom plus external field" is prepared 
in the eigenstate |a) of the Hamiltonian ifo = H &tora + iZfi e id- Then interaction H mt between 
atomic and photon subsystems is switched on. Our aim is to describe possible outcomes of 
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this event. The partial rates of the decay of the initial states \a) into various open channels 
\/3) are given by the expressions 



P 



2ir\T? a (E)\ 2 dp(E), (8) 



where p(E) denotes the densit y o f states in the final state, and the transition operator T 
satisfies the operator equation |3J| 

f(E) = F int + F int (l - P a )— \-f (E), (9) 

E — H 

Here P a is a projection operator on the initial state \a). In Eq. (jHJ) both the matrix element 
of the transition operator and the density of states are to be computed at the energy, 
corresponding to the shifted energy of the initial state \a): E = E a + AE a Finally, the 
energy shift AE a of the final state is related to the diagonal matrix element of the transition 
operator via implicit equation 

AE a = ReT aa (E) (10) 

The form of the latter equation suggests a possiblity of an iterative scheme for the determi- 
nation of the level shift. Equations (jH|)- (fTU|) provide the basis of our calculation scheme for 
evaluating the rates of various multiphoton processes. 

In the alternative formalism j^lj], when the MPI process is treated as a scattering phe- 
nomenon, the expression for the T-operator must be corrected for the virtual transitions 
between various photodetachment channels. This correction comes from a careful exami- 
nation of the boundary conditions which must be imposed on a scattering wave function 
in the Floquet representation. In its simplest form, this can be achieved by dividing the 
conventional expression for the matrix element of the T-operator by a Bessel function Jo 
The problem outlined by Eqs. (jHJ) - (fTU|) is a decay, or an initial state, problem. As 
such, it requires an accurate evaluation of the initial state which subsequently decays into 
various open channels. Indeed, Eq. (jHJ) is a result of evaluation of the quantum mechani- 
cal amplitude (f3\^f(t)) where \I/(t) evolves from the state \a) at t — |3j|. Therefore no 
additional correction to the transition amplitudes Q is needed. 
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By introducing in Eq. (JSJ) a complete set of states of Hq, this equation can be rewritten 
in the form of a spectral representation: 

Here the sign of ie gives the rule of bypassing the pole when performing the integration over 
the continuum spectrum. 

The shift of the energy of the initial state is explicitely included into the set of Eqs.(|lljl. 
This circumstance may be especially important for the atomic systems with more that one 
electron. When the escaping electron leaves the atomic system, the remaining ion still 
interacts with the electromagnetic field, the energy shift in the Eqs. flll)) takes into account 
effects of this interaction. 

In Fig. Elwe give a graphical representation of Eq. fTT| . Here a straight line with an 
arrow to the right represents a target atom and the dashed lines are used to depict the 
photons. A vertex with two electron lines and one photon line represents a dipole matrix 
element (|IJ) which incorporates many-electron correlation in the target. For simplicity we 
do not show the quadrupole and monopole matrix elements (fTTj) which are parametrically 
small for not very strong fields F < 1. A rectangular block stands for the T-matrix (jllj) . In 
the low-field regime the integral term in the right-hand side of (jllj) can be ignored and the 
atomic ionization is described by the bare matrix element (a, n\Hi nt \b, n±l). The strong field 
effects are incorporated in the integral term and include multiple absorption and emission 
of the laser photons. 

The sum over the spectrum of Hq in includes summation over various number of 
photons as well as summation over bound atomic states and integration over the continuous 
spectrum of the atom. Computation of the integral in (jllj) can be greatly simplified by 
introducing a discrete set of target pseudostates which provides an adequate quadrature 
rule. As the result of such a discretization the set of equations (fTTj) becomes a linear system 
on the unknown elements of the T-matrix. Once this linear system is solved, all information 
about the integral and differential features of the MPI process can be obtained from the 
matrix elements of the T-operator taken at the energy corresponding to the shifted energy 
of the final state as in Eq. (jHj) (the so-called on-shell matrix elements). We note that in the 
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m=n+l 



FIG. 1: Diagram representation of the Lippmann-Schwinger equation for multiphoton ionization. 
The graphical symbols are described in the text. 

weak field limit summation in (1111) can be restricted to the atomic variables and we arrive 

n 

to the expression for the dipole matrix element given by Eq. (12) of Kheifets and Bray [35]. 



A. Square well model 

To illustrate feasibility of our approach we first apply it to a model problem of an MPI 
process from a one-dimensional square well. It is this model system that was considered 



by Burke et al. |24| in their seminal paper which gave rise to the spectacular success of 
the i?-matrix Floquet theory. We consider here an electron bound initially in a square well 
potential V = —2.5 a.u. for < x < 1 and V — for x > 1, with the boundary condition 
R{0) = imposed on the wave functions. This potential supports only one bound state a 
with an energy E = —0.4657 a.u. We consider an MPI process when the electric field with 
the frequency uj = 0.2 is applied such that at least three photons are needed to ionize the 
system. To solve this problem we follow the steps outlined in Sec. [H] The set of equations 
(II lj) is converted into a linear system on the T-matrix elements by choosing a suitable 
discretization procedure for the continuous spectrum integration. The summation over a set 
of intermediate states 7 in the Eq. is a sum over various numbers of photons n 7 and 
summation and integration over the set of variables specifying the field-free electron states. 
For the square well model the latter are determined by a single quantum number (momentum 
of the states belonging to the continuous spectrum or energy of the discrete spectrum). For a 
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given number of photons in the set of intermediate states 7 we first determine if there is a pole 
in the i nteg ral over the continuous spectrum. Similarly to electron scattering calculations 
of Bray [34], we divide the momentum integration interval (0, 00) into five subintervals. 

If the pole is present in the integral, the first two subintervals are chosen to be (0, fc po i e ) 
and (fcpoie, 2A; po i e ) with a typical number of 60 momentum points in each interval. Then the 



delta-function singularity is iso 
by a modified Simpson's rule 



ated and the remaining principle value integral is evaluated 
The remaining part of the momentum integral is divided as 
follows: (2fcp i e , 4) (50 integration points), (4, 10) (40 points) and (10, 350) (40 points). These 
intervals are pole-free and the integration is performed by using a conventional Simpson's 
rule. If there is no pole in a given channel the integration procedure remains essentially the 
same except for the boundary of the first and second intervals which is chosen at one atomic 
unit of momentum. In this case the conventional Simpson's rule is used throughout. 

We retained various numbers of photons in the intermediate state of Eq. (fTTj) . For 
simplicity, we count this number from the baseline of three photons in the initial state 
which we denote as a — \a, 3). In the intermediate and final states we can then have 
negative numbers of photons. For example, n = — 1 in the final state I7) would mean that 
four photons have been absorbed. Using this convention, the calculations we performed can 
be denoted as (n min , n max ) meaning that the number of photons in the intermediate and final 
states ranges from n mm to n max . 

The problem which arises immediately when one tries to apply Eqs. (111)1 for the square- 
well model is the singularity of the matrix elements when both states a,/3 lie in the 
continuum. This is, of course, a consequence of the choice of the gauge of the vector potential 
in formulas (JHJ)- Fortunately, as we shall show in the section, this problem can be avoided for 
real atomic systems by choosing the interaction Hamiltonian in the Kramers-Hennerberger 
(acceleration) form. For the square-well potential such choice of the interaction Hamiltonian 
is not possible due to the singular character of the potential. We shall have, therefore, to 
devise a suitable procedure of the regularization of the matrix elements of operator @ when 
both states belong to continious spectrum. There is vast amount of literature devoted to 
th. - dy - such matr , x e,e m en ts for ™ ous potentlals (e., We dealt „* „. 
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problem along the lines described in the paper by Mercouris et al. |23j]. One can show that 
for any potential for which the asymptotic behavior of the continuous wave function is given 
by Rk oc sin (kr + 5), the following result holds: 

] R ki (r)^R kf (r) dr = V(h,k f ) + °° S { \~ 6f) P— L. (12) 
o r * / 

Here V(ki, kf) is a regular function and symbol P has a usual meaning of the principal 

value integral. Note that this divergence is present only when the velocity form (as in the 

present paper) or the length form are used for the electromagnetic interaction operator, it is 

absent if the acceleration form is used. Presence of the singularity in the matrix elements in 

the velocity and length gauges is a consequence of the use of the spectral representation for 

the Green function in and should disappear from the final result after the summation 

over the whole spectrum is carried out. A convenient way to deal with such integrals is 

to introduce a suitable regularization procedure. Different procedures of this sort can be 

devised. We shall illustrate the use of the two simplest methods. 

The first regularization scheme consists in computing all the divergent integrals over some 
finite interval (0, D), hoping that final results depend only very weakly upon D if its value 
is suitably chosen. That this is indeed so can be seen from the results for the level shift 
presented in (JTJ). According to Eq. (JTUJ), this quantity is related to the real part of the 
diagonal matrix element of the T-operator. 

To determine the shift we devised a simple iterative procedure, based on the Eq. ()1U|). The 
linear system (fTT|) was solved with a trial value AE a = 0, than the value given by the r.h.s 
of the Eq. (fTUj) was adopted as a new trial value for AE a , the procedure was repeated untill 
convergence was achieved. Typically, three iterations were sufficient to achieve convergence 
on the level of a fraction of a percent. 

The results of a (-4,4) calculation using this strategy are shown in the Table H] The 
calculation was performed for different values of the cutoff parameter D. As one can see, the 
values for the level shift are fairly stable with respect to the variations of the regularization 
parameter D. The stability of the results with D becomes somewhat worse for very strong 
fields (F = 0.2 a.u. in the Table HJ). This, however, is to be expected as for such strong 
fields the (-4,4) calculation is probably too small and a larger number of photons should 
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be included to obtain more accurate results. Results in the Table |U should be compared 



with the plot illustrating the level shift as a function of the e 



ectric field strength from 



42|. Although the numerical 



the time-dependent i?-matrix calculation of Burke and Burke 
values are not reported in this paper, a visual comparison of our data with the the plot quite 
satisfactory. 



TABLE I: Dependence of the level shift from the (-4,4) calculation upon the regularization pa- 
rameter D. 



D 



F = 0.06 a.u. 



F = 0.1 a.u. 



F = 0.14 a.u. 



F = 0.2 a.u. 



30 
40 
50 
60 
70 



-0.023307 
-0.023305 
-0.023311 
-0.023316 
-0.023313 



-0.065753 
-0.065704 
-0.065794 
-0.065743 
-0.065749 



-0.13132 
-0.13045 
-0.13245 
-0.13386 
-0.13193 



-0.28124 
-0.27645 
-0.29102 
-0.29239 
-0.29839 



Another regularization scheme, which we found to give somewhat more stable results for 
the ionization rates is based on the regularization formula for a principal value integral: 

p i^ f =H^ + k (k t -k f y (13) 

Again, one can hope that if the regularization is properly implemented, the final results do 
not depend upon the regularization parameter e. We analyze these results in the form of 
the partial ionization rates computed according to Eq. (JSjl. The momentum kf of the final 
state is determined via the energy of the initial state a, number of photons in the final state 
rif, and the energy shift AE, given by Eq. (JTUJ): 

k 2 

= E a + AE — n f u> (14) 

The total ionization rate is the sum of the partial rates over all the open channels. 

In the Table ITT1 we show the partial ionization rates r^r^Ts from our (-3,3) calculation. 
The data presented in the Table clearly indicate convergence as e — > 0. In the Table 11111 
we present the partial ionization rates Tj where % has a meaning of the number of absorbed 
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TABLE II: Convergence of the partial ionization rates T% T^, from the (-3,3) calculation with 
respect to the regularization parameter e in 1131 The field strength is F = 0.1 a.u. 



e x 10 2 T 3 (10" 4 a.u.) 



100 


1.582 


50 


1.829 


25 


2.045 


12.5 


2.216 


6.25 


2.334 


3.125 


2.404 


1.563 


2.441 


0.781 


2.454 


0.390 


2.456 


0.195 


2.460 



r 4 (10~ 6 a.u.) 


r 5 (lo 


1.7958 


5.584 


0.9439 


5.866 


0.2198 


6.165 


0.1252 


6.518 


0.4667 


6.552 


0.7188 


6.529 


0.8496 


6.547 


0.9187 


6.578 


0.9537 


6.610 


0.9682 


6.622 



photons for the field strength F = 0.1 a.u. We used different number of laser photons 
(n min ,n max ), with n max = 3,4 and n min = -4,-3,-2,-1,0. In these calculations the 
number of open channels may differ. For example, in the (-4,3) calculation one may have 
the final states with three, four, five, six or seven photons absorbed. We remind the reader 
that we count the number of photons from the baseline of three photons in the initial state. 

TABLE III: Ionization rates for F = 0.1 a.u. 

T 3 (10~ 4 a.u.) T 4 (10" 6 a.u.) T 5 (10~ 6 a.u.) 

-2,3 2.4586 1.0628 6.7172 

-3,3 2.4562 0.9537 6.6101 

-4,3 2.4569 0.9644 6.5025 

-2,4 2.3228 0.8676 6.6303 

-3,4 2.3208 0.7817 6.5132 



As one can see, the partial rates corresponding to processes with large number of photons 
absorbed (four and five) show certain stability with respect to the number of photons in- 
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eluded in the calculation, which gives us confidence in the results obtained for these partial 
rates. One may note, in particular, that for F = 0.1 a.u. the partial rate T 5 becomes larger 
than the rate IY 




0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 

F (a.u.) 



FIG. 2: Branching ratios of the partial ionization rate to the total rate for various open channels 
in the square well model. 

The partial contribution of various open channels to the total ionization rate can be 
judged from the data presented in the Fig. El where we plot the branching ratios for the 
channels with various numbers of photons in the final state. Again, one can see that for not 
very large field strengths ( F < 0.05 a.u.) the contributions of the channels with more than 
three photons absorbed (labeled n = 4 and n = 5 in the Figure) can be neglected, though 
for larger fields their effect is becoming increasingly important. 

Incidentally, one may note, that even (0,3) calculation would not correspond exactly to 
the third order perturbation theory. Indeed, solution of the coupled set of equations (fTTj) 
amounts effectively to summation of an infinite subset of the perturbation theory terms as 
is illustrated in the Fig. [T] 

The situation becomes somewhat more complicated for stronger fields (F > 0.1 a.u.) 
where the level shift is so large that the channel closing may occur. A look at the data for 
the shift presented in the Table H] shows that when F « 0.14 a.u. the energy of the initial 
bound state a is: E ~ —0.6 a.u. and absorption of the three photons is no longer sufficient 
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to ionize the system. A more precise value of the field for which the channel corresponding 
;o absorption of the three photons closes is F ~ 0.1435 a.u. as shown by Burke and Burke 



42j. Rich variety of phenomena may occur when field approaches this critical value (and 



larger ones, corresponding to closing of other channels), in particular ionization rates as 



42|. 



functions of electric field intensity may have discontinuities 

In the Fig. El we plot the total ionization rate as a function of the field strength from 
the (-4,3) calculation. The same data for the total ionization rate are presented in the Table 
] In the figure we make a comparison with the time-dependent R-matrix calculation by 



w 



et al. 



rich was found to be almost identical to the earlier Floquet R-matrix results by Burke 
24|. 




FIG. 3: Total ionization rate as a function of the field strength. The solid line shows the present 
(-4,3) calculation. Results of the time-dependent ii-matrix calculation of Burke and Burke 
shown as points. 



42| are 



Our data show certain trend to deviate from those obtained in Burke and Burke 
and Burke et al. 



42] 



24j for smaller field intensities, around F = 0.02 a.u. This deviation 
could probably be attributed to difficulties of numerical character which are due to the 
necessity to manipulate the divergent integrals. These are most difficult to cope with for 
smaller fields when rates are small. This difficulties may easily be avoided for 'real world' 
problems. For any realistic atomic system, writing the electromagnetic interaction operator 
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in the acceleration form would make all the integrals well-defined and convergent. 

TABLE IV: Total ionization rate from (-4,3) calculation. 



F(a.u.) r(a.u.) F(a.u.) r(a.u.) 

0.025 0.192xl0~ 7 0.14 0.237xl0 _1 

0.05 0.135xl0~ 5 0.15 0.119xl0~ 2 

0.075 0.291xl0" 4 0.16 0.358xl0~ 2 

0.1 0.254xl0" 4 0.175 0.356X10" 1 

0.125 0.279x10^2 



B. Hydrogen atom 



We turn now to a more realistic system - the hydrogen atom in the presence of a 
monochromatic linearly-polarized laser field. As we mentioned above, this system is in some 



sense simpler than the model square- wel 
form of the interaction Hamiltonian 



problem, as we can use the Kramers-Hennebeger 



44j, which, as we shall see, allows to avoid all the 



problems of divergence of the matrix elements. 

If electromagnetic field is treated classically, the Kramers-Henneberger Hamiltonian can 
be obtained from the minimal-coupling Hamiltonian H min given by Eqs. (J2HH1) by means of 
a canonical transformation Hkk = €, lT H m i n e~ lT generated by the operator (44]: 



1 t t 
f = —J A(r)p dr + —j A 2 (r) dr 
o o 



(15) 



where A is the vector potential and p is the momentum operator. 

In the present approach, when electromagnetic field is described quantum-mechanically, 
we need a quantum analog of the transformation (fL3|) . Such transformation is known as 



the Pauli-Fierz canonical transformation 



451 ] . For the case of a sufficiently intense linearly 



polarized monochromatic light, when all the processes of spontaneous radiation are neglected 
and the quantized vector potential is given by the Eq.(0J), with onl y th e laser mode photons 
preserved in the expansion, this transformation assumes the form 4q : 
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r = r' + -Q 

P = p' 

Q = Q' 

p = p,_ v (16) 

to 



where a = 2eJ — , e is polarization vector, operators Q' and P' are expressed in terms of 
creation and annihilation operators in Eg. ({%])) via: 



a k = -j=(P + iQ 



-^( p -*Q) ( 17 ) 



k y/2 

In Eqs. fll6|) . (fT7j) primed quantities refer to original (laboratory) frame, and non-primed 
to the Kramers-Henneberger frame. 

Under the transformation (|16J) the minimal-coupling Hamiltonian (J2J becomes: 

p 2 1 

-f^KH — -7. 1" -Afield + H int KK , (18) 

2 r 

where -f/fi e id has the same meaning as in Eq. (0) and the interaction Hamiltonian has the 
form: 

AntKH = - 1 1 ~, + ~ , (19) 

\r + ct\ r 

where the operator a = F joj 1 is related to the operator of the electric field intensity F. 

Certain amount of care should be taken when physically relevant information is to be 
extracted from the transformed Hamiltonian Hkr, in particular, when our aim is to study 
time-evolution of a quantum system prepared at the moment t = in a given eigenstate 
of the field-free Hamiltonian. The Kramers-Henneberger transformation consists, basically, 
in transforming the Schrodinger equation into a non-inertial frame moving with electron 
oscillating in the field (in the quantum version of this transformation this can be seen 
from the first of Eqs. ljlfij) . the quantity on the r.h.s of this equation being essentially an 
electric field operator). Both the initial and final state vectors must then be transformed 
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accordingly 



47 . . This circumstance may introduce additional complexity in the problem. 



Below, we report results on the total ionization rate for multiphoton ionization of hydrogen 
atom. One can present the following argument showing that if we are interested in integral 
characteristics, such as total ionization rate, one need not transform initial and final states. 
Indeed, solution of the set of equations (jllj) with the Kramers- Henneberger Hamiltonian 
(fT%j) . (|Tnj) is essentially a study of the singular points of the resolvent operator (E — Hku) 1 
(we will be looking, of course, for the points lying on the unphysical sheets of the Riehmann 
surface). Total ionization rate from the ground state of hydrogen obtained below from 
such a study of the Hamiltonian Hkk is in fact related to an imaginary part of one such 
singular point. On the other hand, total ionization rate obtained using the minimal-coupling 
Hamiltonian is a singular point of the resolvent operator {E — -f^ m in) 1 which is connected 
to the resolvent operator (E — H^n)^ 1 by means of the canonical transformation (|16|). The 
latter transformation does not change position of singular points, hence, study of singular 
points Hkk alone (by means of solving set of Eqs. lfTT]) ) should give us correct value for the 
total ionization rate. 

To apply the strategy based on Eqs. fjllj) to the system described by the Kramers- 
Hennebeger Hamiltonian, one should be able to compute the matrix elements of the operator 
(jT9|) . The procedure which we devised for this purpose is described below. 



1. Calculation of matrix elements of Hint- 

We shall need the following formula for the matrix elements of the operator of the elec- 
tromagnetic field, which can be obtained analogously to the Eqs.©: 

(a, n\aF\b, m) m — (a\aF\b) n — > oo, \m — n\ = 1 

(a,n\aF\b,m) =0 \m — n\ ^ 1 (20) 

In Eq. (|2Tlj) F is the classical vector of electric field strength, a is an arbitrary vector operator 
not acting on the photon variables. The correction term to the non-zero matrix element in 
Eq. fl20l) is of the order of and is negligible in the limit of large intensities which is of 
interest to us in the present paper. 
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Consider first the matrix element 

M mp =(n + p\(aF) m \n) (21) 

We are interested now only in photonic variables, therefore we omitted in (J21j) any refer- 
ence to atomic variables. It is easy to see that M mp ^ only if m = p + 2k with an integer 
nonnegative k. If this condition is satisfied, it is easy to show using Eqs.((H) and simple 
combinatorial analysis that 

( aF\ m m\ 

M mv = -, r— , r- (22) 

p \ 2 J ( m + p y. ( m 2 P )- 

We can now compute a more complicated matrix element 

(n + p\ exp {aF}\n} (23) 
By expanding the exponential function in the above expression, with the use of Eq. (|21jl and 



known series expansion for the Bessel function I p (x) 49] 



we can obtain the following expression: 

(n + p\exp{aF}\n) = I p (aF), (25) 

where I p (x) is a modified Bessel function of order p [49]. 

We may now turn to the computation of the matrix elements of the operator (J3*|). Rep- 
resenting (JHJ) as a Fourier transform: 

1 If e iq ( r+6 ^ 



, dq, (26) 
a | 2n 2 J q 2 

using expression f!25[) and known integral representation for the modified Bessel function 



I p {x) 



49|: 



I p (x) = - J e xcose cos p6 d6 (27) 
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one can obtain the following formula: 

1 



n + p 



n 



cos p9 



1 



\n + F/oo 

To proceed further, we use a spherical harmonics expansion 

1 



ri + Fcos6/uj 2 \ 



<:W 



(28) 



E 

k=0 



— ; utt \— sienicosi 

2fc + lr^ +lL V 



l fc Y k0 (r 



(29) 



Ti + Fcos6/lu 2 

where r< (r>) is the smaller (greater) of rt and F cosd/u 2 . Here the field F is directed along 
the z-axis. Eq. (|2U|) allows separation of the radial and angular variables. Angular parts are 
evaluated analytically using integrals of products of several spherical functions 



larpc 
U: 



/ 



Y hmi (n)Y l2m2 (n)Y hma (n) do 



(2Z 1 + l)(2Z 2 + l)(2Z 3 + l) fh k h\(h k h\ 



V 4vr VO 

Radial integrals are computed numerically. 



,mi m 2 m 3 



5. Numerical Results 

We first illustrate our procedure of the calculation of the matrix elements of the interaction 

Hamiltonian. With the help of the formulas (J28)) and (J29|) we calculated the integrals (a, n + 

p\H- mt \b,n) . In the case of hydrogen atom the atomic states \a) are characterized by energies 

and angular momenta. We computed integrals for sufficiently dense grid of energies. As 

for the range of momenta and a number of photons retained in the calculation, their upper 

limits depend on the value of the electric field strength. For moderately high field strengths 

(F « 0.1 a.u.) it was sufficient to compute integrals with I < Z max = 3 and p < p max = 7. 

The integrals were computed and stored on a disk. 

As one can see from the Eq. ()28|) , in the weak field limit p = 1 and the interaction Hamil- 

Fr 

tonian in the Kramers-Hennerberger form reduces to the operator — — - which is commonly 
used in the first order calculations performed in the acceleration gauge. The parameter 
which measures the departure of the matrix element of operator (|28j) from the first order 
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result is the ratio F/u 2 . This departure can be quite significant. To give an illustration 
of the relative role of the higher order corrections in Eq. (J28|) we present in Table |3 few 
values of the matrix elements (ls,n = l\Hi n t\2p, n = 0) and (ls,n = l\Hi nt \kp, n = 0) with 
momentum k = 1 a.u. for different frequencies and electric field strength of F = 0.0534 a.u. 
First order matrix elements are also given in the Table IVl (marked as PT). As can be seen, 
for small frequencies the deviation of the matrix elements from the first order values can be 
significant. 



TABLE V: Matrix elements of the operator ©) for F=0.0534 a.u. 



Angular 


(Is, n 


= 11^12^,^ = 0) 


(Is, n = 


1 \H- mt \kp,n = 0) 


frequency 


Present 


PT 


Present 


PT 


0.65 


0.03070 


0.03118 


0.02126 


0.02143 


0.184 


0.04678 


0.08261 


0.11335 


0.26796 



To facilitate comparison with the literature, we performed a series of calculations for the 
same set of laser field strengths and frequencies as reported by Dorr et al. who used a 
combination of the Floquet ansatz and the i?-matrix methods. Our results for the level shift 
and total ionization rate together with the literature values are summarized in the Table fVTl 

When presenting their data for the level shifts, Dorr et al. |27j who performed their 
calculation in the length gauge subtracted the ponderomotive energy Ep = F 2 /4uj 2 from 
the total shift. Since the level shifts in the acceleration and the length gauge are related as 
AE A = AE L - E P Q the two sets of data are directly comparable. 

We consider first the case of the laser field with u = 0.65 a.u. and F = 0.0534 a.u. 
As noted by Dorr et al. (^J, this is still a perturbative regime. In Eqs. (lllJ) we need not, 
therefore, retain very large set of intermediate states 7 to achieve convergence. 

In the sum over the intermediate states we included the hydrogen atom states with 
orbital momentum 1 — 0,1 and a set of photonic states with numbers of photons ranging 
from n min = —2 to n max = 3. Counting the photons we use here the same convention we 
used in the previous section, when dealing with the square well potential model. We count 
the total number of photons in the system from the minimal number of photons needed to 
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TABLE VI: Total ionization rates and shifts for atomic hydrogen in a linearly polarized field of 
frequency ui and strength F (a.u.). Numbers in parenthesis represent data of Dorr et al. [^tJ 



Angular 


Field 


Total ionization rate 




Shift 


frequency 


strength 


Present 


Dorr et al. [27J 


Present 


Dorr et al. [27J 


0.65 


0.0534 


0.00263 


0.00256 


0.000364 


0.000360 




0.754 


0.13 


0.14 


0.176 


0.195 


0.184 


0.0169 


9.2 x 10~ 6 


8.8 x 10~ 6 


-0.002910 


-0.002543 




0.0534 


1.33 x 10~ 3 


1.40 x 10~ 3 


-0.0280 


-0.0257 



ionize the initial state. Thus, for the angular frequency uo = 0.65 and the ground state of 
the hydrogen atom, we have one photon in the initial state. The state with n = — 1 will 
then correspond to the state in which two photons have been absorbed. 

As one can see from the Table IVI| we achieve quite a good agreement with the data of 
Dorr et al. [^tJ, both for the total rate and level shift. For much larger field and the same 
frequency u = 0.65 a.u. we need a larger set of intermediate states to be included Eq. (fTTf . 
Convergence was achieved when we included the hydrogen atom states with I < / max = 4 
and photonic states with numbers of photons ranging from n m j n = — 2 to n ma , x = 3. 

For the process of genuine multiphoton ionization (a; = 0.184) we also have reasonably 
good agreement with the results of the Floquet R-matrix calculation [^J. For the field 
strength F = 0.0169 we used atomic states with angular momenta I < 4 and photonic states 
with n m i n = — 1 and n max = 5. This field value is still in the perturbative regime. 

Truly nonperturbative effects appear for larger field strengths. To obtain the values of 
the shift and total rate for F = 0.0534 we used atomic states with I < 4 and photonic states 
with n m - m = —2 and n max = 5. 



III. CONCLUSION 

We developed a non-perturbative formalism which describes the MPI process by a set of 
coupled integral equations of the Lippmann-Schwinger type. In contrast to the scattering 
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formalism given in the literature 



17 



19(| , we do not rely on the field modified Volkov 



states. Instead, we employ a more convenient set of the field-free atomic states which is 
particularly advantageous for complex atomic systems with more than one target electron. 
This approach, however, can only be realized if one is able to generate a complete set of 
discrete target pseudostates providing an accurate quadrature rule. In this respect we rely 
on the CCC method which demonstrated its ability to build an efficient pseudostate basis 
for one and two electron targets. This gives us confidence that we can implement our 
computational scheme for nonperturbative description of MPI in complex atomic systems 
such as the helium atom. Thus we should be able to extend a very successful application of 
the CCC method to the weak-filed double photoionization to the strong field domain. 

We demonstrated utility of our approach for a model problem of a square well potential 
and a hydrogen atom. Also, we demonstrated that the partial ionization rates, computation 
of which is usually the most difficult part of the problem, remain stable with respect to the 
number of photons retained in the calculation. 
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